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Abstract —We consider a compressive hyperspectral imaging 
reconstruction problem, where three-dimensional spatio-spectral 
information about a scene is sensed by a coded aperture snapshot 
spectral imager (CASSI). The CASSI imaging process can be 
modeled as suppressing three-dimensional coded and shifted vox¬ 
els and projecting these onto a two-dimensional plane, such that 
the number of acquired measurements is greatly reduced. On the 
other hand, because the measurements are highly compressive, 
the reconstruction process becomes challenging. We previously 
proposed a compressive imaging reconstruction algorithm that is 
applied to two-dimensional images based on the approximate 
message passing (AMP) framework. AMP is an iterative al¬ 
gorithm that can be used in signal and image reconstruction 
by performing denoising at each iteration. We employed an 
adaptive Wiener filter as the image denoiser, and called our 
algorithm “AMP-Wiener.” In this paper, we extend AMP-Wiener 
to three-dimensional hyperspectral image reconstruction, and 
call it “AMP-3D-Wiener.” Applying the AMP framework to the 
CASSI system is challenging, because the matrix that models the 
CASSI system is highly sparse, and such a matrix is not suitable 
to AMP and makes it difficult for AMP to converge. Therefore, we 
modify the adaptive Wiener filter and employ a technique called 
damping to solve for the divergence issue of AMP. Our approach 
is applied in nature, and the numerical experiments show that 
AMP-3D-Wiener outperforms existing widely-used algorithms 
such as gradient projection for sparse reconstruction (GPSR) 
and two-step iterative shrinkage/thresholding (TwIST) given a 
similar amount of runtime. Moreover, in contrast to GPSR and 
TwIST, AMP-3D-Wiener need not tune any parameters, which 
simplifies the reconstruction process. 

Index Terms —Approximate message passing, CASSI, com¬ 
pressive hyperspectral imaging, gradient projection for sparse 
reconstruction, image denoising, two-step iterative shrink- 
age/thresholdng, Wiener filtering. 


1. Introduction 

A. Motivation 

A hyperspectral image is a three-dimensional (3D) image 
cube comprised of a collection of two-dimensional (2D) 
images (slices), where each 2D image is captured at a specific 
wavelength. Hyperspectral images allow us to analyze spectral 
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information about each spatial point in a scene, and thus 
can help us identify different materials that appear in the 
scene [1]. Therefore, hyperspectral imaging has applications 
to areas such as medical imaging [2,3], remote sensing [4], 
geology [5], and astronomy [6]. 

Conventional spectral imagers include whisk broom scan¬ 
ners, push broom scanners [7,8], and spectrometers [9]. In 
whisk broom scanners, a mirror refiects light onto a single 
detector, so that one pixel of data is collected at a time; 
in push broom scanners, an image cube is captured with 
one focal plane array (FPA) measurement per spatial line of 
the scene; and in spectrometers, a set of optical bandpass 
filters are tuned in steps in order to scan the scene. The 
disadvantages of these techniques are that (i) data acquisition 
takes a long time, because they require scanning a number of 
zones linearly in proportion to the desired spatial and spectral 
resolution; and (n) large amounts of data are acquired and 
must be stored and transmitted. For example, for a megapixel 
camera (10^ pixels) that captures a few hundred spectral 
bands (> 100 spectral channels) at 8 or 16 bits per frame, 
conventional spectral imagers demand roughly 10 megabytes 
per raw spectral image, and thus require space on the order of 
gigabytes for transmission or storage, which exceeds existing 
streaming capabilities. 

To address the limitations of conventional spectral imaging 
techniques, many spectral imager sampling schemes based on 
compressive sensing [10-12] have been proposed [13-15]. The 
coded aperture snapshot spectral imager (CASSI) [13,16-18] 
is a popular compressive spectral imager and acquires image 
data from different wavelengths simultaneously. In CASSI, the 
voxels of a scene are first coded by an aperture, then dispersed 
by a dispersive element, and finally detected by a 2D FPA. 
That is, a 3D image cube is suppressed and measured by a 
2D array, and thus CASSI acquires far fewer measurements 
than those acquired by conventional spectral imagers, which 
significantly accelerates the imaging process. In particular, for 
a data cube with spatial resolution of A x M and L spectral 
bands, conventional spectral imagers collect MNL measure¬ 
ments. In contrast, CASSI collects measurements on the order 
of M{N + T — 1). Therefore, the acquisition time, storage 
space, and required bandwidth for transmission in CASSI are 
reduced. On the other hand, because the measurements from 
CASSI are highly compressive, reconstructing 3D image cubes 
from CASSI measurements becomes challenging. Moreover, 
because of the massive size of 3D image data, it is desirable 
to develop fast reconstruction algorithms in order to realize 
real time acquisition and processing. 
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Fortunately, it is possible to reconstruct the 3D cube from 
the 2D measurements according to the theory of compres¬ 
sive sensing [10-12], because the 2D images from different 
wavelengths are highly correlated, and the 3D image cube 
is sparse in an appropriate transform domain, meaning that 
only a small portion of the transform coefficients have large 
values. Approximate message passing (AMP) [19] has recently 
become a popular algorithm that solves compressive sensing 
problems, owing to its promising performance and efficiency. 
Therefore, we are motivated to investigate how to apply AMP 
to the CASSI system. 

B. Related work 

Several algorithms have been proposed to reconstruct image 
cubes from measurements acquired by CASSI. First, the recon¬ 
struction problem for the CASSI system can be solved by ^i- 
minimization. In Arguello and Arce [20], gradient projection 
for sparse reconstruction (GPSR) [21] is utilized to solve for 
the -minimization problem, where the sparsifying transform 
is the Kronecker product of a 2D wavelet transform and a 
ID discrete cosine transform (DCT). Besides using ^i-norm 
as the regularizer, total variation is a popular alternative; 
Wagadarikar et al. [16] employed total variation [22,23] as 
the regularizer in the two-step iterative shrinkage/thresholding 
(TwIST) framework [24], a modified and fast version of 
standard iterative shrinkage/thresholding. Apart from using the 
wavelet-DCT basis, one can sparsify image cubes by dictio¬ 
nary learning [14], or using Gaussian mixture models [25]. An 
interesting idea to improve the reconstruction quality of the 
dictionary learning based approach is to use a standard image 
with red, green, and blue (RGB) components of the same 
scene as side information [14]. That is, a coupled dictionary 
is learned from the joint datasets of the CASSI measurements 
and the corresponding RGB image. We note in passing that 
using color sensitive RGB detectors directly as the FPA of 
CASSI is another way to improve the sensing of spectral 
images, because spatio-spectral coding can be attained in a 
single snapshot without requiring extra optical elements [26]. 

Despite the good results attained with the algorithms men¬ 
tioned above, they all need manual tuning of some param¬ 
eters, which may be time consuming. In GPSR and TwIST, 
the optimal regularization parameter could be different in 
reconstructing different image cubes. In dictionary learning 
methods, although the parameters can be learned automatically 
by methods such as Markov Chain Monte Carlo, the learning 
process is usually time consuming. Moreover, the patch size 
and the number of dictionary atoms in dictionary learning 
methods must be chosen carefully. 

C. Contributions 

In this paper, we develop a robust and fast reconstruction 
algorithm for the CASSI system using approximate message 
passing (AMP) [19]. AMP is an iterative algorithm that can ap¬ 
ply image denoising at each iteration. Previously, we proposed 
a 2D compressive imaging reconstruction algorithm, AMP- 
Wiener [27], where an adaptive Wiener filter was applied as 
the image denoiser within AMP. Our numerical results showed 


that AMP-Wiener outperformed the prior art in terms of both 
reconstruction quality and runtime. The current paper extends 
AMP-Wiener to reconstruct 3D hyperspectral images from 
the CASSI system, and we call the new approach ‘AMP-3D- 
Wiener.” Because the matrix that models the CASSI system is 
highly sparse, structured, and ill-conditioned, applying AMP 
to the CASSI system becomes challenging. For example, (i) 
the noisy image cube that is obtained at each AMP iteration 
contains non-Gaussian noise; and (n) AMP encounters diver¬ 
gence problems, i.e., the reconstruction error may increase 
with more iterations. Although it is favorable to use a high- 
quality denoiser within AMP, so that the reconstruction error 
may decrease faster as the number of iteration increases, we 
have found that in such an ill-conditioned imaging system, 
applying aggressive denoisers within AMP causes divergence 
problems. Therefore, besides using standard techniques such 
as damping [28,29] to encourage the convergence of AMP, 
we modify the adaptive Wiener filter and make it robust to 
the ill-conditioned system model. There are existing denoisers 
that may outperform the modified adaptive Wiener filter in a 
single step denoising problem. However, the modified adaptive 
Wiener filter fits into the AMP framework and allows AMP 
to improve over successive iterations. 

Our approach is applied in nature, and the convergence of 
AMP-3D-Wiener is tested numerically. We simulate AMP- 
3D-Wiener on several settings where complementary random 
coded apertures (see details in Section IV-A) are employed. 
The numerical results show that AMP-3D-Wiener reconstructs 
3D image cubes with less runtime and higher quality than other 
compressive hyperspectral imaging reconstruction algorithms 
such as GPSR [21] and TwIST [16,24] (Figure 3), even 
when the regularization parameters in GPSR and TwIST have 
already been tuned. These favorable results provide AMP-3D- 
Wiener major advantages over GPSR and TwIST. First, when 
the bottleneck is the time required to run the reconstruction 
algorithm, AMP-3D-Wiener can provide the same reconstruc¬ 
tion quality in 100 seconds that the other algorithms provide 
in 450 seconds (Figure 3). Second, when the bottleneck is 
the time required for signal acquisition by CASSI hardware, 
the improved reconstruction quality could allow to reduce the 
number of shots taken by CASSI by as much as a factor 
of 2 (Figure 8). Finally, the reconstructed image cube can 
be obtained by running AMP-3D-Wiener only once, because 
AMP-3D-Wiener does not need to tune any parameters. In 
contrast, the regularization parameters in GPSR and TwIST 
need to be tuned carefully, because the optimal values of these 
parameters may vary for different test image cubes. In order 
to tune the parameters for each test image cube, we run GPSR 
and TwIST many times with different parameter values, and 
then select the ones that provide the best results. 

The remainder of the paper is arranged as follows. We 
review CASSI in Section II, and describe our AMP based 
compressive hyperspectral imaging reconstruction algorithm 
in Section III. Numerical results are presented in Section IV, 
while Section V concludes. 
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11. Coded Aperture Snapshot Spectral Imager 
(CASSI) 

A. Mathematical representation of CASSI 

The coded aperture snapshot spectral imager (CASSI) [18] 
is a compressive spectral imaging system that collects far 
fewer measurements than traditional spectrometers. In CASSI, 
(i) the 2D spatial information of a scene is coded by an 
aperture, (n) the coded spatial projections are spectrally shifted 
by a dispersive element, and (m) the coded and shifted pro¬ 
jections are detected by a 2D FPA. That is, in each coordinate 
of the FPA, the received projection is an integration of the 
coded and shifted voxels over all spectral bands at the same 
spatial coordinate. More specifically, let /o(x, A) denote the 
voxel intensity of a scene at spatial coordinate {x^y) and at 
wavelength A, and let T{x,y) denote the coded aperture. The 
coded density T{x,y)fo{x^y,X) is then spectrally shifted by 
the dispersive element along one of the spatial dimensions. The 
energy received by the FPA at coordinate (x, y) is therefore 

g{x,y)= [ T{x,y - S{X))fo{x,y - S{X),X)dX, (1) 
Ja 

where S{X) is the dispersion function induced by the prism at 
wavelength A. Suppose we take a scene of spatial dimension 
M by and spectral dimension L, i.e., the dimension of the 
image cube is M x A" x L, and the dispersion is along the 
second spatial dimension y, then the number of measurements 
captured by the FPA will be + — 1). If we approximate 
the integral in (1) by a discrete summation and vectorize the 
3D image cube and the 2D measurements, then we obtain a 
matrix-vector form of (1), 

g = Hfo + z, (2) 

where fo is the vectorized 3D image cube of dimension 
n = MNL, vectors g and z are the measurements and the 
additive noise, respectively, and the matrix H is an equivalent 
linear operator that models the integral in (1). In this paper, we 
assume that the additive noise z is independent and identically 
distributed (i.i.d.) Gaussian. With a single shot of CASSI, the 
number of measurements is m = M{N L — 1), whereas K 
shots will yield m = KM{N L — 1) measurements. The 
matrix H in (2) accounts for the effects of the coded aperture 
and the dispersive element. A sketch of this matrix is depicted 
in Figure 1(a) when K = 2 shots are used. It consists of a set 
of diagonal patterns that repeat in the horizontal direction, each 
time with a unit downward shift, as many times as the number 
of spectral bands. Each diagonal pattern is the coded aperture 
itself after being column-wise vectorized. Just below, the next 
set of diagonal patterns are determined by the coded aperture 
pattern used in the subsequent shot. The matrix H will thus 
have as many sets of diagonal patterns as FPA measurements. 
Although H is sparse and highly structured, the restricted 
isometry property [30] still holds, as shown by Arguello and 
Arce [31]. 

B. Higher order CASSI 

Recently, Arguello et al. [32] proposed a higher order model 
to characterize the CASSI system with greater precision, and 



(a) The matrix H for standard CASSI 



(b) The matrix H for higher order CASSI 

Figure 1: The matrix H is presented for iC = 2,M = A = 8, 
and L — A. The circled diagonal patterns that repeat horizontally 
correspond to the coded aperture pattern used in the first FPA shot The 
second coded aperture pattern determines the next set of diagonals. In 
(a) standard CASSI, each FPA shot captures M(A-hL — 1) = 88 
measurements; in (b) higher order CASSI, each FPA shot captures 
M{N -hL-hl) = 104 measurements. 


improved the quality of the reconstructed 3D image cubes. In 
the standard CASSI system model, each cubic voxel in the 3D 
cube contributes to exactly one measurement in the FPA. In 
the higher order CASSI model, however, each cubic voxel is 
shifted to an oblique voxel because of the continuous nature of 
the dispersion, and therefore the oblique voxel contributes to 
more than one measurement in the FPA. As a result, the matrix 
H in (2) will have multiple diagonals as shown in Figure 1(b), 
where there are sets of 3 diagonals for each FPA shot, account¬ 
ing for the voxel energy impinging into the neighboring FPA 
pixels. In this case, the number of measurements with K = 1 
shot of CASSI will be m = M{N + L -f- 1), because each 
diagonal entails the use of M more pixels (we refer readers 
to [32] for details). 

In Section IV, we will provide promising image recon¬ 
struction results for this higher order CASSI system. Using 
the standard CASSI model, our proposed algorithm produces 
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similar advantageous results over other competing algorithms. 
III. Proposed Algorithm 

The goal of our proposed algorithm is to reconstruct the 
image cube fo from its compressive measurements g, where 
the matrix H is known. In this section, we describe our 
algorithm in detail. The algorithm employs (i) approximate 
message passing (AMP) [19], an iterative algorithm for com¬ 
pressive sensing problems, and (n) adaptive Wiener filtering, a 
hyperspectral image denoiser that can be applied within each 
iteration of AMP. 


divergence [35] between the estimated and the true posteriors 
subject to a first order and a second order moment matching 
constraints between fo and Hfo [36]. If the measurement 
matrix H is i.i.d. Gaussian and the empirical distribution of 
fo converges to some distribution on M, then the sequence 
of the mean square error achieved by AMP at each iteration 
converges to the information theoretical minimum mean square 
error asymptotically [37]. 

Moreover, if the matrix is i.i.d. random, then the noise in 
the scalar channel (3) can be viewed as asymptotically i.i.d. 
Gaussian [19,34,38]. 


A. Image denoising in scalar channels 

Below we describe that the linear imaging system model 
in (2) can be converted to a 3D image denoising problem 
in scalar channels. Therefore, we begin by defining scalar 
channels, where the noisy observations q of the image cube fo 
obey 

q = fo + V, (3) 

and V is the additive noise vector. Recovering fo from q is 
known as a 3D image denoising problem. 


B. Approximate message passing 

Algorithm framework: AMP [19] has recently become a 
popular algorithm for solving signal reconstruction problems 
in linear systems as defined in (2). The AMP algorithm 
proceeds iteratively according to 


f*+i =r?t(H^r* + f*), (4) 

r* = g - Hf* + + f*"!)), (5) 

where is the transpose of H, R = m/n represents the 
measurement rate, r]t{') is a denoising function at the t-th 
iteration, 7][{s) = ^r?t(s), and (u) = for 

vector u = {ui^U2,... ^Un). We will explain in Section III-E 
how f^ and are initialized. The last term in (5) is called 
the “Onsager reaction term” [19,33] in statistical physics. 
This Onsager reaction term helps improve the phase transition 
(trade-off between the measurement rate and signal sparsity) of 
the reconstruction process over existing iterative thresholding 
algorithms [19]. In the t-th iteration, we obtain the estimated 
image cube f^ and the residual r^. We highlight that the 
vector -f f^ in (4) can be regarded as a noise-corrupted 
version of fo in the t-th iteration with noise variance 
and therefore rit{’) is a 3D image denoising function that is 
performed on a scalar channel as in (3). Let us denote the 
equivalent scalar channel at iteration t by 

q*=H^r*+r = fo+v*, (6) 


where the noise level is estimated by [34], 



i=l 


(7) 


and r\ denotes the i-th component of the vector in (5). 

Theoretical properties: AMP can be interpreted as min¬ 
imizing a Gaussian approximation of the Kullback-Leibler 


C. Damping 

We have discussed in Section III-B that many mathematical 
properties of AMP hold for the setting where the measurement 
matrix is i.i.d. Gaussian. When the measurement matrix is not 
i.i.d. Gaussian, such as the highly structured matrix H defined 
in (2), AMP may encounter divergence issues. A standard 
technique called “damping” [28,29] is frequently employed 
to solve for the divergence problems of AMP, because it only 
increases the runtime modestly. 

Specifically, damping is an extra step within AMP iterations. 
In (4), instead of updating the value of f^+^ by the output of 
the denoiser + f^), we assign a weighted average 

of 77t(H^r^ + f^) and f^ to f^+^ as follows, 

=a-r/t(H^r* + f*) + (l-a)-f*, (8) 

for some constant 0 < a < 1. Similarly, after obtaining 
in (5), we add an extra damping step that updates the value 
of to be (T • -h (1 — (t) • where the value of a is the 
same as that in (8). 

AMP has been proved [28] to converge with sufficient 
damping, under the assumption that the prior of fo is i.i.d. 
Gaussian with fixed means and variances throughout all iter¬ 
ations, and the amount of damping depends on the condition 
number of the matrix H. Note that other AMP variants [29,39, 
40] have also been proposed in order to encourage convergence 
for a broader class of measurement matrices. 

D. Adaptive Wiener filter 

We are now ready to describe our 3D image denoiser, which 
is the function in the first step of AMP iterations in (4). 

Sparsifying transform: Recall that in 2D image denoising 
problems, a 2D wavelet transform is often performed, and 
some shrinkage function is applied to the wavelet coefficients 
in order to suppress noise [41,42]. The wavelet transform 
based image denoising method is effective, because natural 
images are usually sparse in the wavelet transform domain, 
i.e., there are only a few large wavelet coefficients and the 
rest of the coefficients are small. Therefore, large wavelet 
coefficients are likely to contain information about the image, 
whereas small coefficients are usually comprised mostly of 
noise, and so it is effective to denoise by shrinking the small 
coefficients toward zero and suppressing the large coefficients 
according to the noise variance. Similarly, in hyperspectral 
image denoising, we want to find a sparsifying transform such 
that hyperspectral images have only a few large coefficients 
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in this transform domain. Inspired by Arguello and Arce [20], 
we apply a wavelet transform to each of the 2D images 
in a 3D cube, and then apply a discrete cosine transform 
(DCT) along the spectral dimension, because the 2D slices 
from different wavelengths are highly correlated. That is, the 
sparsifying transform ^ can be expressed as a Kronecker 
product of a DCT transform ^ and a 2D wavelet transform 
W, i.e., ^ ^ (g) W, and it can be shown that ^ is an 

orthonormal transform. Let denote the coefficients of 
in this transform domain, i.e., Our 3D image 

denoising procedure will be applied to the coefficients 
Besides 2D wavelet transform and ID DCT, it is also possible 
to sparsify 3D image cubes by dictionary learning [14] or 
Gaussian mixture models [25]. Moreover, using an endmember 
mixing matrix [43] is an alternative to DCT for characterizing 
the spectral correlation of 3D image cubes. In this work, 
we focus on a 2D wavelet transform and ID DCT as the 
sparsifying transform, because it is an efficient transform that 
does not depend on any particular types of image cubes, and an 
orthonormal transform that is suitable for the AMP framework. 

Parameter estimation in the Wiener filter: In our previous 
work [27] on compressive imaging reconstruction problems 
for 2D images, one of the image denoisers we employed 
was an adaptive Wiener filter in the wavelet domain, where 
the variance of each wavelet coefficient was estimated from 
its neighboring coefficients within a 5 x 5 window, i.e., the 
variance was estimated locally. 

As an initial attempt, we applied the previously proposed 
AMP-Wiener to the reconstruction problem in the CAS SI 
system defined in (2). More specifically, the previously pro¬ 
posed adaptive Wiener filter is applied to the noisy coefficients 
6>q. Unfortunately, AMP-Wiener encounters divergence issues 
for the CAS SI system even with significant damping such 
SiS a = 0.01 in (8). AMP-Wiener diverges, because it is 
designed for the setting where the measurement matrix is 
i.i.d. Gaussian, whereas the measurement matrix H defined 
in (2) is highly structured and not i.i.d., and we found in 
our numerical experiments that the scalar channel noise 
in (6) is not i.i.d. Gaussian. On the other hand, because the 
Wiener filter allows to conveniently calculate the Onsager term 
in (5), we are motivated to keep the Wiener filter strategy, 
although the scalar channel (6) does not contain i.i.d. Gaussian 
noise. Seeing that estimating the coefficient variance from its 
neighboring coefficients (a 3 x 3 or 5 x 5 neighboring window) 
does not produce reasonable reconstruction for the CAS SI 
system, we modify the local variance estimation to a global 
estimation within each wavelet subband. The coefficients 0^ 
of the estimated (denoised) image cube f^ are obtained by 
Wiener filtering, which can be interpreted as the conditional 
expectation of Of given 6>q under the assumption of Gaussian 
prior and Gaussian noise. 


max|0,Pfv. 




iPl-a!)+a! 


(^q,i - Pi,t) + 

(^q,i ~ Mi,t) + Mi,t, (9) 


where 6>q ^ is the i-th element of 0^, and and are 


the empirical mean and variance of 6>q • within an appropriate 
wavelet subband, respectively. Taking the maximum between 
0 and ensures that if the empirical variance 

of the noisy coefficients is smaller than the noise variance , 
then the corresponding noisy coefficients are set to 0. After ob¬ 
taining the denoised coefficients 0\, the estimated image cube 
in the t-th iteration satisfies f^ = ^~^0\ = ^^0\. Therefore, 
the adaptive Wiener filter as a denoiser function can be 
written as 

f‘+i = r?*(q*) 

= (max{0, Vi - a?I}Vi-i (^q* - + tit) , 

( 10 ) 

where 0 is a zero matrix, Vt is a diagonal matrix with z/?^ 
on its diagonal, I is the identify matrix, is a vector that 
contains and max{', •} is operating entry-wise. 

We apply this modified adaptive Wiener filter within AMP, 
and call the algorithm “AMP-3D-Wiener.” We will show in 
Section IV that only a moderate amount of damping is needed 
for AMP-3D-Wiener to converge. 


E. Derivative of adaptive Wiener filter 

The adaptive Wiener filter described in Section III-D is 
applied in (4) as the 3D image denoising function 
The following step in (5) requires p'P'), i.e., the derivative 
of ?7t(-). We now show how to obtain pp’). It has been 
discussed [27] that when the sparsifying transform is or¬ 
thonormal, the derivative calculated in the transform domain 
is equivalent to the derivative in the image domain. According 
to (9), the derivative of the Wiener filter in the transform 
domain with respect to is max{0, Because 

the sparsifying transform ^ is orthonormal, the Onsager term 
in (5) can be calculated efficiently as 


1 


M(q‘)) = ;E 

iGX 




( 11 ) 




where X is the index set of all image cube elements, and the 
cardinality of X is n = MNL. 

We focus on image denoising in an orthonormal transform 
domain and apply Wiener filtering to suppress noise, because 
it is convenient to obtain the Onsager correction term in (5). 
On the other hand, other denoisers that are not wavelet- 
DCT based can also be applied within the AMP framework. 
Metzler et al. [44], for example, proposed to utilize a block 
matching and 3D filtering denoising scheme (BM3D) [45] 
within AMP for 2D compressive imaging reconstruction, and 
run Monte Carlo [46] to approximate the Onsager correction 
term. However, the Monte Carlo technique is accurate only 
when the scalar channel (6) is Gaussian. In the CAS SI system 
model (2), BM4D [47] may be an option for the 3D image 
denoising procedure. However, because the matrix H is ill- 
conditioned, the scalar channel (6) that is produced by AMP 
iterations (4,5) is not Gaussian, and thus the Monte Carlo 
technique fails to approximate the Onsager correction term. 

Having completed the description of AMP-3D-Wiener, we 
summarize AMP-3D-Wiener in Algorithm 1, where fAMP 
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denotes the image cube reconstructed by AMP-3D-Wiener. 
Note that in the first iteration of Algorithm 1, initialization 
of and may not be necessary, because is an all-zero 
vector, and the Onsager term is 0 at iteration 1. 


Algorithm 1 AMP-3D-Wiener 


Inputs: g, H, a, maxiter 
Outputs: fAMP 

Initialization: f^ = 0, = 0 

for t = 1 : maxiter do 

1) r*=g-Hf‘ + ir*-iiEr=i 

^t-1 


max{0,zy2^_i-a2_ j 


2) = a • + (1 — a) • 

3) -|- f^ 

4) = 

5) 6>^ = ^q* 

7) f*+i = a ■ ^'^01 + (1 - a) • P 

end for 

fAMP = 


IV. Numerical Results 


In this section, we provide numerical results where we 
compare the reconstruction quality and runtime of AMP- 
3D-Wiener, gradient projection for sparse reconstruction 
(GPSR) [21], and two-step iterative shrinkage/thresholding 
(TwIST) [16,24]. In all experiments, we use the same coded 
aperture pattern for AMP-3D-Wiener, GPSR, and TwIST. In 
order to quantify the reconstruction quality of each algorithm, 
the peak signal to noise ratio (PSNR) of each 2D slice in the 
reconstructed cubes is measured. The PSNR is defined as the 
ratio between the maximum squared value of the ground truth 
image cube fo and the mean square error of the estimation f, 
i.e.. 


PSNR = 10 • log 


10 


in.^^x,y,X 




\Ea:,y,A (/(a: ,y,X) fo,(x,y,X) 


where f(x,y,x) denotes the element in the cube f at spatial 
coordinate (x, y) and spectral coordinate A. 

In AMP, the damping parameter a is set to be 0.2. Recall 
that increasing the amount of damping helps prevent the 
divergence of AMP-3D-Wiener, and that the divergence issue 
can be identified by evaluating the values of from (7). 
We select 0.2 as the damping parameter value, because 0.2 
is the maximum damping value such that AMP-3D-Wiener 
converges in all the image cubes we test. The divergence 
issues of AMP-3D-Wiener can be detected by evaluating 
the value of obtained by (7) as a function of iteration 
number t. Recall that estimates the amount of noise in 
the noisy image cube q^ at iteration t. If AMP-3D-Wiener 
converges, then we expect the value of to decrease as t 
increases. Otherwise, we know that AMP-3D-Wiener diverges. 
The choice of damping mainly depends on the structure of 
the imaging model in (2) but not on the characteristics of the 


image cubes, and thus the value of the damping parameter a 
need not be tuned in our experiments. 

To reconstruct the image cube fo, GPSR and TwIST mini¬ 
mize objective functions of the form 

?=arginini||g-Hf ||2 + /3 ■</>(£), (12) 

where 0(-) is a regularization function that characterizes the 
structure of the image cube fo, and (3 is sl regularization 
parameter that balances the weights of the two terms in the 
objective function. In GPSR, 0(f) = ||^f||i; in TwIST, the 
total variation regularizer is employed, 

L M N . 

A=1 x=l y=l 

\ 1/2 

+{f{x,y + l,\)-f{x,y,X)f\ . (13) 

Note that the role of the ^i-norm of the sparsify- 
ing coefficients in GPSR is to impose the overall spar¬ 
sity of the sparsifying coefficients, whereas the to¬ 
tal variation in TwIST encourages spatial smoothness in 
the reconstructed image cubes. The implementation of 

GPSR is downloaded from “http://www.lx.it.pt/ mtf/GPSR/,” 
and the implementation of TwIST is downloaded from 
‘ ‘http ://w ww. disp. duke, edu/proj ects/C AS Sl/experimentaldata/ 
index.ptml.” The value of the regularization parameter (3 
in (12) greatly affects the reconstruction results of GPSR and 
TwIST, and must be tuned carefully. We select the optimal 
values of 0 for GPSR and TwIST manually, i.e., we run 
GPSR and TwIST with 5 — 10 different values of (3, and select 
the results with the highest PSNR.^ The typical value of the 
regularization parameter for GPSR is between 10“^ and 10“^, 
and the value for TwIST is around 0.1. We note in passing that 
the ground truth image cube is not known in practice, and 
estimating the PSNR obtained using different 0 may be quite 
involved and require oracle-like information when using GPSR 
and TwIST. Reweighted -minimization [48] does not need 
regularization parameter tuning, and has been shown to out¬ 
perform -minimization by Candes et al. [48]. However, the 
existing reweighted -minimization implementations require 
either QR decomposition [49] of the measurement matrix H 
or the null space of H, which requires H to be expressed 
as a matrix. That said, H is a very large matrix, and we 
implement it as a linear operator. Therefore, implementing the 
reweighted -minimization that is applicable to the system 
model in (2) is beyond the scope of this paper, and the 
reweighted -minimization is not included in our simulation 
results. There exist other hyperspectral image reconstruction 
algorithms based on dictionary learning [14,25]. In order to 
learn a dictionary that represents a 3D image, the image cube 
needs to be divided into small patches, and the measurement 
matrix H also needs to be divided accordingly. Dividing the 

^As an example, we simulate GPSR with many different values for /3, and 
obtain that for /3 = 1 • 10“^, 2 • 10“^, 3 • 10“^, 4 • 10“^, 5 • 10“^, 6 • 10“^, 
and 7 • 10“^, the corresponding PSNRs of the reconstructed cubes are 31.25 
dB, 32.30 dB, 32.82 dB, 32.99 dB, 33.02 dB, 33.09 dB, and 33.06 dB. 
Therefore, we select /3 = 6 • 10“^ for this specific image cube. We follow 
the same procedure to select the optimal (3 values for each test image cube. 
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Figure 2: The Lego scene. (The target object presented in the experimental results was not endorsed by the trademark owners and it is used here 
as fair use to illustrate the quality of reconstruction of compressive spectral image measurements. LEGO is a trademark of the LEGO Group, 
which does not sponsor, authorize or endorse the images in this paper. The LEGO Group. All Rights Reserved, http://aboutus.iego.com/en- 
us/legal-notice/fair-play/.) 


measurement matrix into smaller patches is convenient for the 
standard CASSI model (Figure 1(a)), because there is a one-to- 
one correspondence between the measurement matrix and the 
image cube, i.e., each measurement is a linear combination 
of only one voxel in each spectral band. In higher order 
CASSI, however, each measurement is a linear combination 
of multiple voxels in each spectral band. Therefore, it is not 
straightforward to modify these dictionary learning methods 
to the higher order CASSI model described in Section II-B, 
and we do not include these algorithms in the comparison. 

A. Test on ''Lego” image cube 

The first set of simulations is performed for the scene shown 
in Figure 2. This data cube was acquired using a wide-band 
Xenon lamp as the illumination source, modulated by a visible 
monochromator spanning the spectral range between 448 nm 
and 664 nm, and each spectral band has 9 nm width. The 
image intensity was captured using a grayscale CCD camera, 
with pixel size 9.9 /im, and 8 bits of intensity levels. The 
resulting test data cube has MxA^ = 256x256 pixels of 
spatial resolution and L = 24 spectral bands. 


Setting 1: The measurements g are captured with K = 
2 shots. The coded aperture in the first shot is generated 
randomly with 50% of the aperture being opaque, and the 
coded aperture in the second shot is the complement of 
the aperture in the first shot. The measurement rate with 
two shots is m/n = KM{N -f L + 1)/{MNL) ^ 0.09. 
Moreover, we add Gaussian noise with zero mean to the 
measurements. The signal to noise ratio (SNR) is defined as 
10 log]^o(/^^/^noise) [20], where jig is the mean value of the 
measurements Hfo and (Jnoise is the standard deviation of the 
additive noise z. In Setting 1, we add measurement noise such 
that the SNR is 20 dB. 

We note in passing that the complementary random coded 
apertures are binary, and can be implemented through pho¬ 
tomask technology or emulated by a digital micromirror 
device (DMD). Therefore, the complementary random coded 
apertures are feasible in practice [20]. Moreover, the comple¬ 
mentary random coded apertures ensure that in the matrix H 
in (2), the norm of each column is similar, which is suitable 
for the AMP framework. However, it is a limitation of the cur¬ 
rent AMP-3D-Wiener that the complementary random coded 
apertures must be employed, otherwise, AMP-3D-Wiener may 
























Figure 3: Runtime versus average PSNR comparison of AMP-3D- 
Wiener, GPSR, and TwIST for the Lego image cube. Cube size is M = 
N = 256, and L = 24. The measurements are captured with K — 2 
shots using complementary random coded apertures, and the number 
of measurements is m = 143,872. Random Gaussian noise is added 
to the measurements such that the SNR is 20 dB. 


diverge. 

Let us now evaluate the numerical results for Setting 1. 
Figure 3 compares the reconstruction quality of AMP-3D- 
Wiener, GPSR, and TwIST within a certain amount of runtime. 
Runtime is measured on a Dell OPTIPLEX 9010 running 
an Intel(R) CoreTM i7-860 with 16GB RAM, and the en¬ 
vironment is Matlab R2013a. In Figure 3, the horizontal axis 
represents runtime in seconds, and the vertical axis is the av¬ 
eraged PSNR over the 24 spectral bands. Although the PSNR 
of AMP-3D-Wiener oscillates during the first few iterations, 
which may be because the matrix H is ill-conditioned, it 
becomes stable after 50 seconds and reaches a higher level 
when compared to the PSNRs of GPSR and TwIST at 50 
seconds. After 450 seconds, the average PSNR of the cube 
reconstructed by AMP-3D-Wiener (solid curve with triangle 
markers) is 26.16 dB, while the average PSNRs of GPSR (dash 
curve with circle markers) and TwIST (dash-dotted curve with 
cross markers) are 23.46 dB and 25.10 dB, respectively. Note 
that in 450 seconds, TwIST runs roughly 200 iterations, while 
AMP-3D-Wiener and GPSR run 400 iterations. 

Figure 4 complements Figure 3 by illustrating the PSNR 
of each 2D slice in the reconstructed cube separately. It is 
shown that the cube reconstructed by AMP-3D-Wiener has 
2 — 4 dB higher PSNR than the cubes reconstructed by GPSR 
and 0.4 — 3 dB higher than those of TwIST for all 24 slices. 

In Figure 5, we plot the 2D slices at wavelengths 488 nm, 
533 nm, and 578 nm in the actual image cubes reconstructed 
by AMP-3D-Wiener, GPSR, and TwIST. The images in these 
four rows are slices from the ground truth image cube fo, the 
cubes reconstructed by AMP-3D-Wiener, GPSR, and TwIST, 
respectively. The images in columns 1 — 3 show the upper- 
left part of the scene, whereas images in columns 4 — 6 
show the upper-right part of the scene. All images are of 
size 128 X 128. It is clear from Figure 5 that the 2D slices 
reconstructed by AMP-3D-Wiener have better visual quality; 
the slices reconstructed by GPSR have blurry edges, and the 



Figure 4: Spectral band versus PSNR comparison of AMP-3D- 
Wiener, GPSR, and TwIST for the Lego image cube. Cube size is 
M — N — 256, and L — 24. The measurements are captured with 
K — 2 shots using complementary random coded apertures, and the 
number of measurements is m — 143,872. Random Gaussian noise is 
added to the measurements such that the SNR is 20 dB. 


slices reconstructed by TwIST lack details, because the total 
variation regularization tends to constrain the images to be 
piecewise constant. 

Furthermore, a spectral signature plot analyzes how the 
pixel values change along the spectral dimension at a fixed 
spatial location, and we present such spectral signature plots 
for the image cubes reconstructed by AMP-3D-Wiener, GPSR, 
and TwIST in Figure 6. Three spatial locations are selected 
as shown in Figure 6(a), and the spectral signature plots 
for locations B, C, and D are shown in Figures 6(b)-6(d), 
respectively. It can be seen that the spectral signatures of 
the cube reconstructed by AMP-3D-Wiener closely resemble 
those of the ground truth image cube (dotted curve with square 
markers), whereas there are discrepancies between the spectral 
signatures of the cube reconstructed by GPSR or TwIST and 
those of the ground truth cube. 

According to the runtime experiment from Setting 1, we 
run AMP-3D-Wiener with 400 iterations, GPSR with 400 
iterations, and TwIST with 200 iterations for the rest of the 
simulations, so that all algorithms complete within a similar 
amount of time. 

Setting 2: In this experiment, we add measurement noise 
such that the SNR varies from 15 dB to 40 dB, which is the 
same setting as in Arguello and Arce [20], and the result is 
shown in Figure 7. Again, AMP-3D-Wiener achieves more 
than 2 dB higher PSNR than GPSR, and about 1 dB higher 
PSNR than TwIST, overall. 

Setting 3: In Settings 1 and 2, the measurements are 
captured with K = 2 shots. We now test our algorithm on 
the setting where the number of shots varies from K = 2 to 
K = 12 with pairwise complementary random coded aper¬ 
tures. Specifically, we randomly generate the coded aperture 
for the k-th shot fork = 1, 3, 5, 7, 9,11, and the coded aperture 
in the {k-f l)-th shot is the complement of the aperture in the 
k-th shot. In this setting, a moderate amount of noise (20 dB) 
is added to the measurements. Figure 8 presents the PSNR of 
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Figure 5: 2D slices at wavelengths 488 nm, 533 nm, and 578 nm in the image cubes reconstructed by AMP-3D-Wiener, GPSR, and TwIST 
for the Lego image cube. Cube size is M — N — 256, and L — 24. The measurements are captured with K — 2 shots using complementary 
random coded apertures, and the number of measurements is m = 143,872. Random Gaussian noise is added to the measurements such that the 
SNR is 20 dB. First row: ground truth; second row: the reconstruction result by AMP-3D-Wiener; third row: the reconstruction result by GPSR; 
last row: the reconstruction result by TwIST. Columns 1 — 3; upper-left part of the scene of size 128 x 128; columns 4 — 6; upper-right part of 
the scene of size 128 x 128. 


the reconstructed cubes as a function of the number of shots, 
and AMP-3D-Wiener consistently beats GPSR and TwIST. 

B. Test on natural scenes 

Besides the Lego image cube, we have also tested our 
algorithm on image cubes of natural scenes [50].^ There are 
two datasets, “natural scenes 2002” and “natural scenes 2004,” 
each one with 8 image data cubes. The cubes in the first 
dataset have L = 31 spectral bands with spatial resolution 
of around 700 x 700, whereas the cubes in the second dataset 
have L = 33 spectral bands with spatial resolution of around 
1000 X1000. To satisfy the dyadic constraint of the 2D wavelet, 
we crop their spatial resolution tobeM = A^ = 512. Because 
the spatial dimensions of the cubes “scene 6” and “scene?” in 
the first dataset are smaller than 512x512, we do not include 
results for these two cubes. 

^The cubes are downloaded from http://personalpages.manchester.ac.uk/staff/ 
d.h.foster/Hyperspectral_images_of_natural_scenes_04.html and http://per- 
sonal pages.manchester.ac.uk/staff/d.h.foster/Hyperspectral_images_of_natural_ 
scenes_02.html. 


The measurements are captured with K = 2 shots, and the 
measurement rate is m/n = KM{N -F L F 1)/{MNL) « 
0.069 for “natural scene 2002” and 0.065 for “natural scene 
2004.” We test for measurement noise levels such that the 
SNRs are 15 dB and 20 dB. The typical runtimes for AMP 
with 400 iterations, GPSR with 400 iterations, and TwIST with 
200 iterations are approximately 2,800 seconds. The average 
PSNR over all spectral bands for each reconstructed cube is 
shown in Tables I and II. We highlight the highest PSNR 
among AMP-3D-Wiener, GPSR, and TwIST using bold fonts. 
It can be seen from Tables I and II that AMP-3D-Wiener usu¬ 
ally outperforms GPSR by 2 — 5 dB in terms of the PSNR, and 
outperforms TwIST by 0.2 — 4 dB, while TwIST outperforms 
GPSR by up to 3 dB for most of the scenes. Additionally, the 
results of 6 selected image cubes are displayed in Figure 9 
in the form of 2D RGB images.^ The four rows of images 
correspond to ground truth, results by AMP-3D-Wiener, results 

^We refer to the tutorial from http://personalpages.manchester.ac.uk/staff/ 
david.foster/Tutorial_HSI2RGB/Tutorial_HSI2RGB.html and convert 3D im¬ 
age cubes to 2D RGB images. 
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Figure 6: Comparison of AMP-3D-Wiener, GPSR, and TwIST on reconstruction along the spectral dimension of three spatial pixel locations 
as indicated in (a). The estimated pixel values are illustrated for (b) the pixel B, (c) the pixel C, and (d) the pixel D. 




Number of shots 


Figure 7: Measurement noise versus average PSNR comparison of 
AMP-3D-Wiener, GPSR, and TwIST for the Lego image cube. Cube 
size is M — N — 256, and L — 24. The measurements are captured 
with K — 2 shots using complementary random coded apertures, and 
the number of measurements is m = 143,872. 


Figure 8: Number of shots versus average PSNR comparison of 
AMP-3D-Wiener, GPSR, and TwIST for the Lego image cube. Cube 
size is N — M — 256, and L — 24:. The measurements are cap¬ 
tured using pairwise complementary random coded apertures. Random 
Gaussian noise is added to the measurements such that the SNR is 20 
dB. 


by GPSR, and results by TwIST, respectively. We can see from 
Figure 9 that the test datasets contain both smooth scenes and 
scenes with large gradients, and AMP-3D-Wiener consistently 
reconstructs better than GPSR and TwIST, which suggests that 
AMP-3D-Wiener is adaptive to various types of scenes. 


V. Conclusion 

In this paper, we considered the compressive hyperspec- 
tral imaging reconstruction problem for the coded aperture 
snapshot spectral imager (CASSI) system. Considering that 
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Figure 9: Comparison of selected image cubes reconstructed by AMP-3D-Wiener, GPSR, and TwIST for the datasets '‘natural scene 2002” 
and "natural scene 2004.” The 2D RGB images shown in this figure are converted from their corresponding 3D image cubes. Cube size is 
]\[ = M = 512, and L = 31 for images in columns 1 — 2 or L = 33 for images in columns 3 — 6. Random Gaussian noise is added to the 
measurements such that the SNR is 20 dB. First row: ground truth; second row: the reconstruction result by AMP-3D-Wiener; third row: the 
reconstruction result by GPSR; last row: the reconstruction result by TwIST. 


SNR 

15 dB 

20 dB 

Algorithm 

AMP 

GPSR 

TwIST 

AMP 

GPSR 

TwIST 

Scene 1 

32.69 

28.10 

31.05 

33.29 

28.09 

31.16 

Scene 2 

26.52 

24.32 

26.25 

26.65 

24.40 

26.41 

Scene 3 

32.05 

29.33 

31.21 

32.45 

29.55 

31.54 

Scene 4 

27.57 

25.19 

27.17 

27.76 

25.47 

27.70 

Scene 5 

29.68 

27.09 

29.07 

29.80 

27.29 

29.42 

Scene 8 

28.72 

25.53 

26.24 

29.33 

25.77 

26.46 


TABLE I: Average PSNR comparison of AMP-3D-Wiener, GPSR, 
and TwIST for the dataset "natural scene 2002” downloaded from [50]. 
The spatial dimensions of the cubes are cropped to M = N = 512, 
and each cube has L = 31 spectral bands. The measurements are 
captured with K — 2 shots, and the number of measurements is 
m = 557,056. Random Gaussian noise is added to the measurements 
such that the SNR is 15 or 20 dB. Because the spatial dimensions of 
the cubes "scene 6” and "scene?” in "natural scenes 2002” are smaller 
than 512 x 512, we do not include results for these two cubes. 


the CAS SI system is a great improvement in terms of imag¬ 
ing quality and acquisition speed over conventional spectral 
imaging techniques, it is desirable to further improve CAS SI 
by accelerating the 3D image cube reconstruction process. 
Our proposed AMP-3D-Wiener used an adaptive Wiener filter 
as a 3D image denoiser within the approximate message 
passing (AMP) [19] framework. AMP-3D-Wiener was faster 
than existing image cube reconstruction algorithms, and also 
achieved better reconstruction quality. 


SNR 

15 dB 

20 dB 

Algorithm 

AMP 

GPSR 

TwIST 

AMP 

GPSR 

TwIST 

Scene 1 

30.48 

28.43 

30.17 

30.37 

28.53 

30.31 

Scene 2 

27.34 

24.71 

27.03 

27.81 

24.87 

27.35 

Scene 3 

33.13 

29.38 

31.69 

33.12 

29.44 

31.75 

Scene 4 

32.07 

26.99 

31.69 

32.14 

27.25 

32.08 

Scene 5 

27.44 

24.25 

26.48 

27.83 

24.60 

26.85 

Scene 6 

29.15 

24.99 

25.74 

30.00 

25.53 

26.15 

Scene 7 

36.35 

33.09 

33.59 

37.11 

33.55 

34.05 

Scene 8 

32.12 

28.14 

28.22 

32.93 

28.82 

28.69 


TABLE II: Average PSNR comparison of AMP-3D-Wiener, GPSR, 
and TwIST for the dataset "natural scene 2004” downloaded from [50]. 
The spatial dimensions of the cubes are cropped to M = N = 512, 
and each cube has L — 33 spectral bands. The measurements are 
captured with K — 2 shots, and the number of measurements is 
m = 559,104. Random Gaussian noise is added to the measurements 
such that the SNR is 15 or 20 dB. 


In AMP, the derivative of the image denoiser is required, 
and the adaptive Wiener filter can be expressed in closed 
form using a simple formula, and so its derivative is easy to 
compute. Although the matrix that models the CASSI system 
is ill-conditioned and may cause AMP to diverge, we helped 
AMP converge using damping, and reconstructed 3D image 
cubes successfully. Numerical results showed that AMP-3D- 
Wiener is robust and fast, and outperforms gradient projec¬ 
tion for sparse reconstruction (GPSR) and two-step iterative 































12 


shrinkage/thresholding (TwIST) even when the regularization 
parameters for GPSR and TwIST are optimally tuned. More¬ 
over, a significant advantage over GPSR and TwIST is that 
AMP-3D-Wiener need not tune any parameters, and thus an 
image cube can be reconstructed by running AMP-3D-Wiener 
only once, which is critical in real-world scenarios. In contrast, 
GPSR and TwIST must be run multiple times in order to find 
the optimal regularization parameters. 

Future improvements: In our current AMP-3D-Wiener al¬ 
gorithm for compressive hyperspectral imaging reconstruction, 
we estimated the noise variance of the noisy image cube within 
each AMP iteration using (7). In order to denoise the noisy 
image cube in the sparsifying transform domain, we applied 
the estimated noise variance value to all wavelet subbands. 
The noise variance estimation and 3D image denoising method 
were effective, and helped produce promising reconstruction. 
However, both the noise variance estimation and the 3D image 
denoising method may be sub-optimal, because the noisy 
image cube within each AMP iteration does not contain i.i.d. 
Gaussian noise, and so the coefficients in the different wavelet 
subbands may contain different amounts of noise. On the other 
hand, in the proposed adaptive Wiener filter, the variances 
of the coefficients in the sparsifying transform domain were 
estimated empirically within each wavelet subband, whereas 
it is also possible to apply Wiener filtering via marginal 
likelihood or generalized cross validation [51]. Therefore, it 
is possible that the denoising part of the proposed algorithm 
can be further improved. The study of such denoising methods 
is left for future work. 

In our current AMP-3D-Wiener, the coded apertures must 
be complementary, because complementary coded apertures 
ensure that the norm of each column in the matrix H in (2) 
is similar, otherwise, AMP-3D-Wiener may diverge. Although 
using complementary coded aperture has practical importance, 
it provides more fiexibility in coded aperture design when such 
a complementary constraint can be removed, and the develop¬ 
ment for AMP-based algorithms without such constraints is 
left for future work. 

Finally, besides reconstructing image cubes from compres¬ 
sive hyperspectral imaging systems, it would also be inter¬ 
esting to investigate problems such as target detection [52] 
and unmixing [53] using compressive measurements from 
hyperspectral imaging systems. We leave these problems for 
future work. 
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